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In M Zhu and Rabitz presented a rapidly convergent iter- 
ative algorithm for optimal control of the expectation value 
of a positive definite observable in a pure-state quantum sys- 
tem. In this paper we generalize this algorithm to a quantum 
statistical mechanics setting and show that it is both efficient 
in the mixed-state case and effective in achieving the control 
objective of maximizing the ensemble average of arbitrary ob- 
servables in the cases studied. 



I. INTRODUCTION 

Much work has recently been done on control of 
pure-state quantum systems using the traditional wave- 
function formalism 10-0]. This work is most important; 
however many physical systems, such as systems initially 
in thermal equilibrium or otherwise described by an en- 
semble of states, or systems where dissipative processes 
are significant, can not be treated using this approach. 
Therefore, a development of optimal control for mixed- 
state quantum systems is necessary. In this paper we 
shall focus on generalizing an efficient iterative algorithm 
for quantum control 0] to a quantum statistical mechan- 
ics setting used in previous work H H . This work is 
closely related to recently published, independently de- 
veloped work by Yukiyoshi, Zhu and Rabitz |^ on quan- 
tum optimal control for systems with dissipation. How- 
ever, in our work we do not consider dissipation terms 
since those terms are represented by non-Hermitian op- 
erators resulting in non-unitary evolution of the system. 
Unfortunately, the very accurate numerical implementa- 
tion of the algorithm we propose depends on unitary evo- 
lution, as do the results on kinematical bounds Q] and 
controllability llfl] , which we use to show that the actual 
global maximum is reached by this algorithm. 



II. MATHEMATICAL SETUP 

As in our previous work, we consider a quantum- 
mechanical system whose state space Ti is a. separable 



Hilbert space. Any mixed state of the system can be rep- 
resented by a density operator p(t) (acting on Ti.) with 
eigenvalue decomposition 



p(i)=^«;fc|*fc(t))(*fc(t)|. 



(1) 



where Wk are the eigenvalues, and \^k{t)) the correspond- 
ing normalized eigenstates of p{t), which evolve in time 
according to the time-dependent Schrodinger equation. 
The eigenvalues satisfy 



< Wfc < 1 Vfc and >, ^k 

k 



(2) 



i.e., they can be ordered in a (possibly finite) non- 
increasing sequence 



Wl > W2 > ■ ■ ■ > Wk > 



>0. 



Unless otherwise mentioned, the word state will in the 
following refer to a mixed state represented by a density 
operator p{t). 

The dynamical law for the system is given by the quan- 
tum Liouville equation 



|p(o=4[^,p(t)]. 



(3) 



where H is the (total) Hamiltonian of the system and 
pita) = po defines the initial state of the system (at time 

t = to). 

Observables are represented by Hermitian operators A 
on Ti. and we define their expectation value to be the 
ensemble average 



(Ait)) = Tr (ip(i)) 



(4) 



The set of bounded linear operators A on H forms 
itself a Hilbert space, usually called Liouville space and 
it is convenient to assign to each operator A (on TC) a 
Liouville ket |^)) denoting its representation in Liouville 
space. The dual of \A)) will be denoted by the Liouville 
bra {{A\. The inner product in Liouville space is defined 
by 



{{A I B)) = Tr (its 



(5) 



Thus, an arbitrary mixed state of the system is repre- 
sented by a LiouviUe ket \p{t))) that satisfies 



^^\pit))) = -^c{tMt))) 



(6) 



with some initial condition \p{to))) — \pa))- C is the 
LiouviUe operator defined by the dual correspondence 



C\p{t)))^[H,p{t)]- 



(7) 



The expectation value {A{t)) of the observable A is given 
by the LiouviUe inner product {{A \ p{t))). 



III. CONTROLLING THE DYNAMICS 



If the number M of external control functions 



m^{h{t).h{t),...jM{t)). 



(8) 



acting on the system is finite and the system is control- 
linear then the total Hamiltonian of the system can be 
decomposed as follows: 



M 



H^Ho 



7 , fm{t)Hn 



(9) 



m— 1 



In this case, the corresponding LiouviUe operator also 
decomposes: 



M 
m— 1 



(10) 



The restrictions imposed on the controls depend on the 
particular system studied. However, a reasonable mini- 
mal requirement for the control functions fm{t) is that 
they should be bounded, measurable, real-valued func- 
tions defined on a time interval [to , tp\ that depends on 
the application. 

In the remainder of this paper we shall furthermore 
assume that there is only one control /(t) acting on the 
system, which is sufficient for many applications of laser 
control. However, we would like to point out that it is 
possible to generalize the algorithm to the case where 
there are multiple controls, such as two laser fields with 
perpendicular polarization driving the system. 

Our goal is to maximize the expectation value (ensem- 
ble average) of a given observable, e.g., the population of 
a particular energy level or subspace of quantum states, 
the energy of a molecular bond, etc., at some fixed target 
time t = tp subject to certain constraints. 

More precisely, we define a functional M,pl|l|] 



WU,Pv.A,)^Wi{p,)-W2{I.P..A,)-W^{!), (11) 



whose value at a certain target time tp we would like to 
maximize. Wi is the expectation value of A which we 
wish to maximize at the target time tp, 



W,U)^{A{tp))^{{A\p,{tp)))- 



(12) 



W2 and W3 are constraint functionals, which we define 
as follows: 



W2{f,Pv,A^) 



tp 



{{A.{t)\§-, + ^C{t)\p,{t)))dt, (13) 



W3{f) = ^ f^f{t)dt. 



(14) 



W2 ensures that the quantum LiouviUe equation is satis- 
fied. W3 constrains the fluence, i.e., the total energy of 
the pulse. 

Pu(i) and Ay{t) are variational trial functions that 
must satisfy the boundary conditions 



Pvito) ^ p{to) ^ Po, Ay{tp) ^ A. 



(15) 



For simplicity we shall in the following choose units such 
that h — I and define dt = ■§[■ 

Eqs (|l2|)-([l5[) are the generalization to LiouviUe space 
of the Hilbert space formulation in [nl . The details of the 
connection with this paper will be discussed in appendix 

0- . . 

The solution of this control problem requires finding 

an admissible control f(t) such that W and thus {A{t)) 
will attain its global maximum at time t = tp. 



IV. ALGORITHM 

We start by guessing an initial control /^"•'(t) and de- 
termining an initial jpl, (t))) by solving 



dM^\t))) 



C, + f'^(t)C, |pW(t))) 



S<^)i 



with initial condition \pv (to))) — \po))- 
For n > 1 and fc = 0, 1 we define 



/(".fc)(t) ^ --((4")(t)|/:i|p("-^-)(t))) (16) 



and solve iteratively 

a,i4")(t))) = -*£("-i)(t)i4")(t))) (18) 

^t\p^:^(t))) = -^6-'°'>{t)\p^^^(t))) (19) 

with the boundary conditions 

l4")(t^))) = |A)), U"^(to))) = |po)). 



We observe that f^^'^\t) is real. Hence £("''^) is Hermi- 

tian and the time-evolution of both \Av (t)) and \pv (t)) 
is unitary, i.e., 



4")(i) 



and 



^(") 



(t) 



= II Poll 



for all t G [ta,tp] and any n. Furthermore, 
\\po\\1=Tt(^pIpo)^Tt{pI)<1. 



(20) 



(21) 



This algorithm can be shown to converge quadratically 
and monotonically as does the pure-state version due to 
Zhu and Rabitz. The details of the proof can be found in 
appendix H. However, we have no guarantee that Wi{f) 
indeed assumes its global maximum for this f{t). Ad- 
ditional criteria, such as kinematical bounds and knowl- 
edge about controllability of the system are necessary to 
decide if the control the algorithm produced is indeed 
optimal in the sense of steering the system to a global 
maximum of Wi{f). 



V. NUMERICAL IMPLEMENTATION 

The differential equations arising from this feedback 
algorithm must be solved numerically. While there are 
many methods of integrating differential equations nu- 
merically, we employ a symmetric split operator method 
[E|,E2l . The main advantage of this method is that it pre- 
serves the norm of the operators involved, which is of 
great importance in this problem. 

We divide the time interval [io,^_F] in subintervals 
[tj, tj-i-i] of a fixed length At — tj+i — tj. On each subin- 
terval [ij,ij+i] we approximate /("''^^(t) by the constant 

Tj = tj + At/2 = tj+i - At/2. (22) 

With this approximation the propagator can be written 
as 

U("''=)(tj+i,i,) = exp{^tAt{£o + &'''\r,)£,)). (23) 

For arbitrary matrices A and B we have 

^-ia(A+B) ^ ^-i(a/2)A^-iaB^-Ha/2)B _ 

up to second order terms in A and B. Thus ( p3| ) agrees 
to second order with 



g-»4i£og-»At/(".'=)(T,)£ig-*Al£o^ 



(24) 



This symmetric splitting is numerically favorable since 
it allows us to reduce the matrix exponentials to a simple 
linear combination of complex exponentials: 



Uo = exp{-iAtCo/2) 

N 

= 5:|a))e-^*/2((a| 

a=l 

ul"^'^(r,)^exp(-zAt/("''=)(T,)/:i) 



N 



Ei^))^ 



-,:At/'"''=)(T,)6 



m 



(25) 



(26) 



b=l 



where \a)) and \b)) are the eigenkets of £o and Ci, respec- 
tively; a and b are the corresponding (real) eigenvalues. 
This leads to 



N 



U(".fc)(r^.) = Y, l((« I &))pe-*'^*'"+''^ ("^■»|a))((a|. 



,,b=l 



(27) 



U^"'''^(t,) agrees up to second order with u'"'''^(ij_i, t^). 
Since £o and £i do not depend on /("■'=), the eigenvalue 
decomposition needs to be done only once, i.e., the only 
quantities that need to be computed in each step of the 
iteration are the complex exponentials g-^^tio^+bf "■ (tj)) 
for all possible values of a and b. 

In order to compute /{tj), we note that 



fit±At)^fit)±Atf^{t) 
to 1st order, and hence we have 



(28) 



At 



+^{(Air){t,)\[C„,£i]pi"\t,))} (29) 



— ((4")(t,)|[/:o,/:i]pi"-^^(t,))). (30) 



VI. ILLUSTRATIVE COMPUTATIONS 

As an example for molecular quantum control, we con- 
sider a Morse oscillator model for a diatomic molecule 
with N discrete energy levels _E„ corresponding to inde- 
pendent vibrational eigenstates \n) of the system. The 
unperturbed Hamiltonian is thus 



N 



i^o = E^'»^ 



(31) 



The interaction Hamiltonian of the driven system can be 
approximated by Hi = f{t)V where f(t) is an external 
laser field that serves as control function, and V is the 
transition operator, which we choose to be of the dipole 
form 



V 



N-1 

E 

n=l 



d„(|n}(n + l| + |n + l)(n|). 



(32) 



This system is completely controllable, which can easily 
be verified using an algorithm described in |10|| . Thus, 
the global minima and maxima of any observable are de- 
termined by the kinematical bounds and these extrema 
are dynamically attainable. 

For the sake of illustration we choose A^ = 4. The cor- 
responding energy levels are Ei — 0.4843, E2 — 1.4214, 
E3 = 2.3691 and E^ = 3.2434 in units of huj^ where 
cjo = 7.8 X VS"^ s-i for HF. 

Let us first assume that the system is initially in the 
ground state, i.e., po = |1)(1| and that our goal is to 
maximize the vibrational energy of the bond, i.e., A — 
Ho- In this case, the results on kinematical bounds in H] 
give 



1.4214 < {A{t)) < 3.2434. 



(33) 



The lower bound is attained exactly if the population of 
level 1 (ground state) is 1. The upper bound is attained 
exactly if the population of level 4 (highest state) is 1. 
Figs 1-3 show the results of our computations using the 
algorithm described above. Starting with a randomly 
generated function / of sufhciently small magnitude and 
A = 4, the observable rapidly approaches its converged 
value within only a few iterations. Fig. 1 shows the final 
pulse f{t), Fig. 2 the corresponding evolution of the pop- 
ulations of energy levels 1 through 4, and Fig. 3 shows 
the evolution of the expectation value of the observable. 
At the target time tp = 200 fs, we observe a nearly com- 
plete inversion of the populations, with the population of 
level four being close to 97%. (A(ij?)} is about 98% of 
the theoretical maximum. 

Secondly, we assume that the system is initially in ther- 
mal equilibrium, i.e., 

N 
n=l 

with weights 

w„ = Cexp{-En/{E4-Ei)). 
This is a Bolzmann distribution with kT = E4 — Ei. 

(J ^ l^-EJkT ^ g-B2/feT ^ ^-Es/kT ^ ^-Ei/kTyl 

is the normalization constant. Concretely, wi = 0.3850, 
W2 = 0.2758, W3 = 0.1976 and w^ = 0.1416. According 
toH, 



the highest population, the second most energetic state 
has the second highest population, etc. Figs 4-6 show the 
results of our computations using the algorithm described 
above. Again, we started with a randomly generated 
function / of sufficiently small magnitude and A = 4. 
Fig. 4 shows the final pulse f{t), Fig. 5 the corresponding 
evolution of the populations of energy levels 1 through 4, 
and Fig. 6 shows the evolution of the expectation value 
of the observable. At the target time tp = 200 fs we 
observe a nearly complete inversion of the populations 
with {A{tF)) being 99% of the theoretical maximum. 
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0.15 




100 150 200 

Time (fs) 
FIG. 1. Optimal pulse for a four- level Morse oscillator with 

Po = |l>(l| 

Populations 
1 




100 150 200 

Time ffsl 
FIG. 2. Evolution of the populations for a four-level Morse 

oscillator with po = |1)(1| 



1.5059 < (A) < 2.2592. 



(34) 



The lower bound is attained in thermal equilibrium. The 
upper bound is attained exactly if the populations are 
inverted, i.e., the most energetic state (here n = 4) has 
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Time (fs) 
FIG. 3. Evolution of the vibrational energy for a four-level 

Morse oscillator with po = |1)(1| 
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Time (fsl 
FIG. 6. Evolution of the vibrational energy for a four-level 

Morse oscillator with po = X^^^j \^)(jA 
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FIG. 4. Optimal pulse for a four-level Morse oscillator with 
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100 150 200 

Time (fsi 
FIG. 5. Evolution of the populations for a four-level Morse 

oscillator with po = ^ _, |^){"-| 



VII. CONCLUSION 

In this paper we demonstrated that an efficient algo- 
rithm for optimal control of quantum systems can be ap- 
plied in a quantum statistical mechanics setting and that 
this algorithm is also highly effective in realizing the con- 
trol objective of maximizing the ensemble average of an 
observable. 
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APPENDIX A: RELATION TO WORK OF 
RABITZ ET AL. 

Our variational functional and Euler-Lagrange equa- 
tions are equivalent to the ones used in |l|] in the pure 
state Hmit, i.e., if /5.„(t) = \^y[t)){ipy{t)\ where \'4'v{t)) is 
a normalized state then 

W = {MtF)\A\MtF)) -ao [^ f{t)dt 

Jto 

-25R / ^XvitWt+iHifMMmt (Al) 

Choose a (time-dependent) complete orthonormal set 
{\ii^n{t)) : n = 1, 2, . . .} such that \-ipi{t)) = \ipv{t)) for aU 
t. Then we have 



Wi = i:i[AMtF) 

n 

= (^,(iF)|Al^«(if )>. 



Furthermore, setting \xv{t)) — Ay{t)\i/jv{t)) we obtain 
{{A,{t) I dtp.it))) = Tr (^A.{t)dtpAt)) 

= J2{Mt)\A,it){dt\Mm{Mt) I Mt)) 

+ Y,{Mt)\Mt)\Mt)){dt{Mm\Mt)) 
= (V'.(t)|i.(t)9t|V'.W) 
+ Y,idt{Mt)\)\Mt)){Mt)\Mt)\Mt)) 

n 

= (v.(t)|A„(i)at|v.(t)) + {dt{Mt)\)Mt)\Mt)) 
= (V.(t)|i„(i)9t|7A„(t)) + {{Mt)\Ar,{t)dt\MtW 

^2^{MtMvit)dt\Mt)) 
= 2Jft(x.(i) I dtMt)) 



and 



{{A,{t)\tC{f,t)p.,{t))) 
= iTr(i,(i)[i?(/,i),p„(t)] 

= E„*(V'n(i)|i.Wi?(/,t)|</>.(i))(V'.(i) I V-nW) 
-E„*(V'n(t)|i.(i)|<^.(i))(V'«WIWi)IV'nW) 

= (x.WN^(/,t) 10.(0) - {iPAt)\iH{f,t)\xAt)) 

= (x.WN^(/,t)|0.(t)) + ((x«WNWOI'^.W))' 

= 2di{xvmH{f,t)\cf,4t)). 



Hence, we have 



W. 



tp 



{{A.{t)\dt+iC{t)\p.{t)))dt 



= 2K / {xv{t)\[dt + iH{f,t)Mv{t))dt 

Jto 

in the pure state case. W3 remains essentially the same, 
i.e., we simply set ao — A/2. The equivalence of the 
Euler-Lagrange equations follows. 



APPENDIX B: PROOF OF CONVERGENCE 
PROPERTIES 



\{{A\pi-\t)))\'< 
as well as 

\fit)\' = 
< 
< 



pi'-Ht) ^<\\Ar,, 



--((4"^(t)|Apl"^(i))) 



< 



1 

A2 

1 

A2 

1 

A2 



4")(i) ^- C^pi^Ht) 

4")(i)[.||A||.||p(")(t) 

4"'(i)[-||A|| 



where ||>Ci|| is the usual operator norm. Thus, 



tp — to 



< \\Ak 



2A 



4"^*) 



11^1 II 



for all n, which establishes the claim. 
Lemma: If U(t, io) satisfies 



dt\J{t,to) ^ ~iC{t)\J{t,to) 



then 



|p(t)))-u(t,io) / \JHt'Mm')))dt' 

is a solution of 

dMm^-^m\p{t))) + m))) 

Proof: Using the product rule and 

dt [ \j\t\to)m')))dt' = v\t,to)m))). 



to differentiate \p{t))) gives 



dt\pit)))^[-iCit)V{t,to)] / V\t',to)m')))dt' 



-\jit,to)vHt,to)mt))) 
-^c{t)\pm + m))). 



After the nth iteration step, the objective functional is 



Theorem 1: Convergence. The sequence {VF^"^} con- 
verges monotonically and quadratically in the control, 
i.e., 



{{A\pi-\tp))) 



[/(-°)(i)]2rfi (Bl) - hm M/("+i) - W(») 



since W^""^ = Ty2(/^"'°\pi"\ ^1"^) = according to Eqs 
® and (|l|). 

Lemma: M^("' is uniformly bounded. 

Proof: Cauchy-Schwarz's inequality and Eqs (20), 
(P) give 



= lim - /''[(5/("+i)(i)]2 + [^/("+i-")(t)]2di (B2) 
Proof: Setting 

\spirHm = \p^:'^'\t)))~\pi"Hm, (bs) 



A r*' 



[/("+l'0)(i)]2-[/("^0)]2di. 



(B4) 



During the iteration 

dtlpi'-Hm = -^[£o + /("■°) WA]|pl"Ht))>- (B5) 
Hence, setting 



and noting that 



(B6) 

(B7) 



/:i/("+i'i)|<5p(")(t))) 

£^j(n+l.l)|^(n+l)(i))) _ /:,/(«+14)|^(")(i))) 

+/:i/("+i'")|p("+i)(0)) - /:i/("+i'i)|pi"+i)(i))) 



we obtain 



5,|<5p(")(t)))=-z/:("+^^i)|5p(")(t))) 

(B8) 



Setting 



exp. 



z / /:("+i^i)(r)d7 



where exp_|_ denotes the time-ordered exponential, the 
formal solution of B8 is (according to the previous 
lemma) given by 



(BIO) 
Observing that 

i4") m = u(i, to, /("+^'i))ut(if , to, /("+^-^))iA)) 

and thus 

((A(")(t)| = ((A|U(tf,to,/^"+^'^^)Ut(t,to,/("+^^^)), 
we arrive at 



{{A\Spi'^\tF))) 
= -I r"((A|U(tf ,to,/^"+'^'^)Ut(t,to,/("+'^^))x 

A|(<^/("+'Vi"+'^ + '5/("+i'"Vi"')(i)))di- 

=-* r((4"+i)(t)ix 

A|(<5/("+^Vi"+'^ + <5/("+i'"Vi"^)(0))di- 

= /*"-z<5/("+i)(t)((4"+i)(t) |/:ip("+i)(t)))dt 
-z<5/("+i'")(t)((4"+i)(t) I /:ipi")(t)))dt 

= A / 5/("+i) (t)/("+i^") (t) + ^/(«+i^") (i)/(«+i^i) (i)dt 

■/to 

= A /*I/("+i^°ni)]' - /("+i^^)(t)/("+^'°ni) 

Jto 

(Bll) 

J^(n+l,n) ^ ^ /"*f/(»+1.0)(i)]2+2/("+l'l)(t)/("+l^°)(t) 

2 Jto 

_^2[/("+i,i)(i)]2 _ 2/("^°)(t)/("+i'i)(t) 

+ [/("'") (t)]2dt 



j^y(n+l)(^)]2_^j^y(n+l,«)(^)]2^^ 



(B12) 



and thus the total variation from n = to up is 

71^ — 1 



(B9) ^T/tA("i=-'0) ^ ^(np) _ ^(0) ^ ^ ^M/("+i' 



n) 



n=0 



"^-'a -■* 



n=0 ^ -^'o 

(B13) 

Since Vl^(") is uniformly bounded, W^"^"^ - VT^^^ is also 
uniformly bounded for all np and thus the sequence 
{(5P^(»-F'.o) . „^ g ]sj^| jg uniformly bounded. 



[^^(«+l)(^)]2 ^ [5/(»+l,«)(^)]2^^ > Q 



for any n implies furthermore that SW^'^'^'^' is an increas- 
ing sequence. Hence, 

lim (5^^"*''°) exists and is finite. 
Consequently 

lim - /''[(5/("+i)(t)]2 + [,5/("+i^")(t)]2rft = 0. 
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